/*******************************************************************************

Name: Figure 2 - persistence analysis

Authors: Jonathan Baker, Lori Bennear, Sheila Olmstead

Input files:
	Primary panel: natlCWS_yr_Panel_mod-small_v2-Match.dta
	alternative panel: natlCWS_yr_Panel_mod.dta

		
*******************************************************************************/


* ========================== START START START ==========================
clear
clear matrix
set more off

* ______________________________________________________________________________
* set paths and create log file
global root "C:/Users/jbaker/Documents/Research/SDWIS/Analysis/STATA"
global input "$root/Input"
global intermediate "$root/Intermediate"
global log "$root/LogFiles"
global output "$root/Output"
global outpanel "$root/OutputPanel"
global figures "$root/Figures"

log using "$log/Figure2_persistence.log", replace name(analysis)

* ============================================================================== 

* Persistence analysis regresssions

* ==============================================================================

// the following block of code borrows from tradeoffAnalysis_v3.do to merge back
//   into the matched dataset the specific violation types. This allows for a
//   definition of microbial violations
	use "$outpanel/natlCWS_yr_Panel_mod.dta", clear
		keep pwsidNUM year HLTH_*
		drop HLTH_999 // matched panel already has HLTH_999
		sort pwsidNUM year
		isid pwsidNUM year
		save "$intermediate/temp_cwspanel.dta", replace
		
	use "$outpanel/natlCWS_yr_Panel_mod-small_v2-Match.dta", clear
		merge 1:1 pwsidNUM year using "$intermediate/temp_cwspanel.dta"
		sort pwsidNUM year
		isid pwsidNUM year
		save "$intermediate/temp_tradeoffPanel.dta", replace

	erase "$intermediate/temp_cwspanel.dta"

	gen microb = HLTH_110 + HLTH_121 + HLTH_122 + HLTH_123 + HLTH_130 + HLTH_140
		// define microbial violations (see tradeoffAnalysis_v3.do)

local id = 2001
local depvar "all microb"
foreach y of local depvar {

	// the following if statement sets the dependent variable
	if "`y'" == "all" {
		local DV = "HLTH_999"
		}
	else if "`y'" == "microb" {
		local DV = "`y'"
		}

display "================= Dependent variable is set to `DV' ================="
// persistence model
local threshold "501 10k 100k"
foreach req of local threshold { 
	// annual dummies (1998 through end of sample): full effect in each year
		reghdfe `DV' T_`req'x1998-T_`req'x`id' ///
		if id`id'==1 & mflag`req'==1 [fweight=freq_`req'], ///
		a(i.pwsidNUM i.state#i.year) vce(cl pwsidNUM)
			display "Model: `e(cmdline)'"
			estimates save "$output/tdummy-sysFEstyrFE-id`id'-`y'-`req'.ster", replace

		reghdfe `DV' T_`req'x1998-T_`req'x`id' Postxsize Postxsize2 ///
		if id`id'==1 & mflag`req'==1 [fweight=freq_`req'], ///
		a(i.pwsidNUM i.state#i.year) vce(cl pwsidNUM)
			display "Model: `e(cmdline)'"
			estimates save "$output/tdummy-sysFEstyrFEfsize2-id`id'-`y'-`req'.ster", replace
}
}		
// Dose-response persistence model: dv = all health violations
reghdfe HLTH_999 ///
	T_501x1998-T_501x2001 ///
	T_10kx1998-T_10kx2001 ///
	T_100kx1998-T_100kx2001 ///
	Postxsize Postxsize2 ///
	if id2001==1 & mflag501==1 [fweight=freq_501], ///
	a(i.pwsidNUM i.state#i.year) vce(cl pwsidNUM)

display "Model: `e(cmdline)'"
estimates save "$output/tdummy-sysFEstyrFEfsize2-id2001-all-DR.ster", replace

// Dose-response persistence model: dv = microbial health violations
reghdfe microb ///
	T_501x1998-T_501x2001 ///
	T_10kx1998-T_10kx2001 ///
	T_100kx1998-T_100kx2001 ///
	Postxsize Postxsize2 ///
	if id2001==1 & mflag501==1 [fweight=freq_501], ///
	a(i.pwsidNUM i.state#i.year) vce(cl pwsidNUM)

display "Model: `e(cmdline)'"
estimates save "$output/tdummy-sysFEstyrFEfsize2-id2001-microb-DR.ster", replace

* ============================================================================== 

* Create Figure 2 plot

* ==============================================================================
set more off
clear all

local depvar "all microb"
foreach dv of local depvar {
	estimates use "$output/tdummy-sysFEstyrFEfsize2-id2001-`dv'-501.ster"
		estimates store `dv'501
	estimates use "$output/tdummy-sysFEstyrFEfsize2-id2001-`dv'-10k.ster"
		estimates store `dv'10k
	estimates use "$output/tdummy-sysFEstyrFEfsize2-id2001-`dv'-100k.ster"
		estimates store `dv'100k
}
* define the psytle for each violation type
local all_p = "p10"
local microb_p = "p3"

* COEFPLOT code
coefplot ///
( ///
all501, keep(T_501x1998) rename(T_501x1998 = "Publish") \ ///
all10k, keep(T_10kx1998) rename(T_10kx1998 = "Mail") \ ///
all100k, keep(T_100kx1998) rename(T_100kx1998 = "Web") label("All Violations") m(O) pstyle(`all_p') ///
) ///
( ///
microb501, keep(T_501x1998) rename(T_501x1998 = "Publish") \ ///
microb10k, keep(T_10kx1998) rename(T_10kx1998 = "Mail") \ ///
microb100k, keep(T_100kx1998) rename(T_100kx1998 = "Web") label("Microbial Violations") m(D) pstyle(`microb_p') ///
), bylabel(1998) || ///
( ///
all501, keep(T_501x1999) rename(T_501x1999 = "Publish") \ ///
all10k, keep(T_10kx1999) rename(T_10kx1999 = "Mail") \ ///
all100k, keep(T_100kx1999) rename(T_100kx1999 = "Web") label("All Violations") m(O) pstyle(`all_p') ///
) ///
( ///
microb501, keep(T_501x1999) rename(T_501x1999 = "Publish") \ ///
microb10k, keep(T_10kx1999) rename(T_10kx1999 = "Mail") \ ///
microb100k, keep(T_100kx1999) rename(T_100kx1999 = "Web") label("Microbial Violations") m(D) pstyle(`microb_p') ///
), bylabel(1999) || ///
( ///
all501, keep(T_501x2000) rename(T_501x2000 = "Publish") \ ///
all10k, keep(T_10kx2000) rename(T_10kx2000 = "Mail") \ ///
all100k, keep(T_100kx2000) rename(T_100kx2000 = "Web") label("All Violations") m(O) pstyle(`all_p') ///
) ///
( ///
microb501, keep(T_501x2000) rename(T_501x2000 = "Publish") \ ///
microb10k, keep(T_10kx2000) rename(T_10kx2000 = "Mail") \ ///
microb100k, keep(T_100kx2000) rename(T_100kx2000 = "Web") label("Microbial Violations") m(D) pstyle(`microb_p') ///
), bylabel(2000) || ///
( ///
all501, keep(T_501x2001) rename(T_501x2001 = "Publish") \ ///
all10k, keep(T_10kx2001) rename(T_10kx2001 = "Mail") \ ///
all100k, keep(T_100kx2001) rename(T_100kx2001 = "Web") label("All Violations") m(O) pstyle(`all_p') ///
) ///
( ///
microb501, keep(T_501x2001) rename(T_501x2001 = "Publish") \ ///
microb10k, keep(T_10kx2001) rename(T_10kx2001 = "Mail") \ ///
microb100k, keep(T_100kx2001) rename(T_100kx2001 = "Web") label("Microbial Violations") m(D) pstyle(`microb_p') ///
), bylabel(2001) || ///
, vertical yline(0) levels(99 95 90) ciopts(recast(rspike rcap rcap))

* EXPORT FIGURE
graph export "$figures/persist_byYear.pdf", replace

* ============================================================================== 

* Create Figure D1 plot

* ==============================================================================-

** annual time dummy plots
set more off
clear all

local depvar "all microb"
foreach dv of local depvar {
	estimates use "$output/tdummy-sysFEstyrFEfsize2-id2001-`dv'-DR.ster"
		estimates store `dv'DR
}
* define the psytle for each violation type
local all_p = "p10"
local microb_p = "p3"

* COEFPLOT code
coefplot ///
( ///
allDR, keep(T_501x1998) rename(T_501x1998 = "Publish") \ ///
allDR, keep(T_10kx1998) rename(T_10kx1998 = "Mail") \ ///
allDR, keep(T_100kx1998) rename(T_100kx1998 = "Web") ///
label("All Violations") m(O) pstyle(`all_p') ///
) ///
( ///
microbDR, keep(T_501x1998) rename(T_501x1998 = "Publish") \ ///
microbDR, keep(T_10kx1998) rename(T_10kx1998 = "Mail") \ ///
microbDR, keep(T_100kx1998) rename(T_100kx1998 = "Web") ///
label("Microbial Violations") m(D) pstyle(`microb_p')), bylabel(1998) || ///
( ///
allDR, keep(T_501x1999) rename(T_501x1999 = "Publish") \ ///
allDR, keep(T_10kx1999) rename(T_10kx1999 = "Mail") \ ///
allDR, keep(T_100kx1999) rename(T_100kx1999 = "Web") ///
label("All Violations") m(O) pstyle(`all_p') ///
) ///
( ///
microbDR, keep(T_501x1999) rename(T_501x1999 = "Publish") \ ///
microbDR, keep(T_10kx1999) rename(T_10kx1999 = "Mail") \ ///
microbDR, keep(T_100kx1999) rename(T_100kx1999 = "Web") ///
label("Microbial Violations") m(D) pstyle(`microb_p')), bylabel(1999) || ///
( ///
allDR, keep(T_501x2000) rename(T_501x2000 = "Publish") \ ///
allDR, keep(T_10kx2000) rename(T_10kx2000 = "Mail") \ ///
allDR, keep(T_100kx2000) rename(T_100kx2000 = "Web") ///
label("All Violations") m(O) pstyle(`all_p') ///
) ///
( ///
microbDR, keep(T_501x2000) rename(T_501x2000 = "Publish") \ ///
microbDR, keep(T_10kx2000) rename(T_10kx2000 = "Mail") \ ///
microbDR, keep(T_100kx2000) rename(T_100kx2000 = "Web") ///
label("Microbial Violations") m(D) pstyle(`microb_p')), bylabel(2000) || ///
( ///
allDR, keep(T_501x2001) rename(T_501x2001 = "Publish") \ ///
allDR, keep(T_10kx2001) rename(T_10kx2001 = "Mail") \ ///
allDR, keep(T_100kx2001) rename(T_100kx2001 = "Web") ///
label("All Violations") m(O) pstyle(`all_p') ///
) ///
( ///
microbDR, keep(T_501x2001) rename(T_501x2001 = "Publish") \ ///
microbDR, keep(T_10kx2001) rename(T_10kx2001 = "Mail") \ ///
microbDR, keep(T_100kx2001) rename(T_100kx2001 = "Web") ///
label("Microbial Violations") m(D) pstyle(`microb_p')), bylabel(2001) || ///
, vertical yline(0) levels(99 95 90) ciopts(recast(rspike rcap rcap))

* EXPORT FIGURE
graph export "$figures/DRpersist_id2001_byYear.pdf", replace
* ========================== DONE DONE DONE ==========================
log close analysis
